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ABSTRACT 

We present a new numerical tool for solving the special relativistic ideal MHD equa¬ 
tions that is based on the combination of the following three key features: (i) a one-step 
ADER discontinuous Galerkin (DG) scheme that allows for an arbitrary order of ac¬ 
curacy in both space and time, (ii) an a posteriori subcell finite volume limiter that is 
activated to avoid spurious oscillations at discontinuities without destroying the nat¬ 
ural subcell resolution capabilities of the DG finite element framework and finally (iii) 
a space-time adaptive mesh refinement (AMR) framework with time-accurate local 
time-stepping. 

The divergence-free character of the magnetic field is instead taken into account 
through the so-called ’’divergence-cleaning” approach. The convergence of the new 
scheme is verified up to 5^^ order in space and time and the results for a set of 
significant numerical tests including shock tube problems, the RMHD rotor and blast 
wave problems, as well as the Orszag-Tang vortex system are shown. We also consider 
a simple case of the relativistic Kelvin-Helmholtz instability with a magnetic field, 
emphasizing the potential of the new method for studying turbulent RMHD fiows. We 
discuss the advantages of our new approach when the equations of relativistic MHD 
need to be solved with high accuracy within various astrophysical systems. 

Key words: magnetohydrodynamics, special relativity, ADER discontinuous Galerkin, 
a posteriori subcell limiter, adaptive mesh refinement, local time stepping 


1 INTRODUCTION 


Special relativistic magnetohydrodynamics (RMHD) is supposed to provide a sufficiently accurate description of the dynamics 
of those astrophysical plasma that move close to the speed of light and which are subject to electromag netic forces dominatin g 
over gravitational f orces. This is the case of high energy astrophysical phenomena li ke extragalac tic jets ( Begelman et al.lll984 b 
gamma-ray bursts ( Kouveliotou et al. 1993) and magnetospheres of neutron stars (|Michellll99lh . In all these physical systems, 
in fact, leaving aside the problem of the origin of relativistic jets, which clearly involves the role of the accretion disc and of 
the corresponding central compact object, general relativistic efiects can be fairly neglected. The degree of complexity related 
to magnetohydrodynamics can of course vary notably. For example, as a first approximation we can neglect dissipation due to 
resistivity, or we can consider the fiuid as a single-component one, although we know for sure that both magnetic reconnection 
and multi-fiuids efiects can become important under specific physical conditions. 

The numerical solution of the special relativistic magnetohydrodynamics equations has been particularly fostered by 
the introduction of Godunov methods based on Riemann solvers, which had alre ady been successfu lly ap plied to r elativis - 
tic hydrodynamics. This was the approach followed in the pioneering works by iKomissarov (|l999l l and iBalsaral (|2001al b 
who implemented for the first time second order Total Variation Diminishing (TVD) schemes with a specific interest to¬ 
wards astrophysical applications. Since then, relativistic magnetohydrodynamics has developed al ong difierent directions 
with impressive results. From one side several approximate Riemann solvers have been introduced ( Mignone Bod3 20061 : 
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Honkkila, h .TanViiiriRn 2007 : Mignone et al. 200d : Kim fc B alsaral 2014 ). From another side, relativistic m a gnetohydrodynam- 


ics has been extended to the general relativistic reg ime 


Baumgarte fc Shapirol 2003 : Anton et al. 20061 : 


Del Zanna et alJ 120071 : iGiacomazzo Rezzollal 120071 ) , and it is currently used to study a variety of high energy physical 


sistive magnetohydrodynamics, with encouraging results (Komissarov 20071: Palenzuela et al. 20091: Dumbser V Zanotti 20091: 

Zenitani et al. 2010 

: Takamoto & Inoud 2011: Bucciantini & Del Zannal 2013lj. Moreover, high order numerical schemes have 

also been pursued ( 

Del Zanna et al.l 

2 OO 3 I: lAnderson et al. 2006lL while simulations of multi-fluids in RMHD are emerging as 

a new frontier (Zenitani et al. l2009l: 

iBarkov et al.l 12014 

h. Finally, Adaptive Mesh Refinement (AMR) within RMHD codes 

has been also considered (jBalsaral 2001bl: Neilsen et al.l 

I 2 OO 6 I: Etienne et al. 2010l: Mignone et al. 2012: Keopens et al. 2012: 


valuable alternative is provided by ADER schemes, which were introduced bv Titarev & Torol ( 

20(4): 

: Toro & Titarev 

(l2nnd 

and became Dooular after the modern reformulation bv Dumbser et al. (|2008l): Dumbser et al. 

( 

|200d 


iBalsara et al. ( 

Q 


In a nutshell, ADER schemes are high order numerical schemes with a sing le step for the time update and they have been 
already applied to the equ ations of relativistic MHD , both in the ideal case (|Dumbser et al.ll2008l : IZanotti fc Dumbserl 120151 ) 
and in the resistive case ( Dumbser fc Zanottil 2009l l . Another common choice that is typically adopted in the majority of 
modern RMHD codes is that of using finite difference or finite volume conservative schemes. Although rather successful, these 
schemes require larger and larger stencils when the order of accuracy is incr eased, a fact that can g ive rise to substantial 
overhead when they are par allelized. Discontinuous Galerkin (DG) schemes ( Gockburn Shu IQSOl : Gockburn et al. IQSOl . 
1990l : Gockburn Shulll998l b on the contrary, do not need any spatial reconstruction and they allow for an arbitrary order of 


accuracy. DG schemes are still relativel y unknown in high energy astrophysics, and only a few investigation s have been per¬ 
formed so far in the relativistic regime ( Zumbusd] 20091 : Radice Rezzoll3 2011 : Zanotti Dumbserl 201 ih . Unfortunately, 
DG schemes suffer from a serious problem, which has negatively affected their popularity. Namely, since they are linear in the 
sense of Godunov’s theorem, they produce oscillations as soon as a discontinuity appears in the solution, even though they 
exploit the conservative formulation of the equations and even though Riemann solvers are used for the computation of the 
fluxes. The procedures that have been adopted to overcome this difficulty can be roughly divided in two classes. From one side, 
it is possible to introduce add i tional numerical diss ipation, either in the form of artificial viscosity Jr. H artmann fc P.HoustonI 
20021 : Persson Peraire 20061 : Gesenek et al. 2013[ ). or by means of filtering ( Radice fc Rezzollal 201 ih . From another side, it 


of nonlinear finite-volume-type slope-limiting procedure (iGockburn & Shu 

19981: Qiu & Shul 20051: 

IBalsara et al.ll2007l: Ij.Zhu et al. 20081: J.Zhu & Oiu 20131: H.Luo et al.l 2007 

: Ih.Krivodonoval 20071). 


either based on nonlinear 

WENO/HWENO reconstruction or by applying a TVB limiter to the higher order moments of the discrete solution. The 
drawback of this strategy is that in most cases the subcell resolution properties of the DG scheme are immediately lost. 


Very recently, a pro mising alternative has been proposed by Dumbser et al. I ll2014lk which is based on a previous idea 


of Glain et al. (2011) and Diot et al. ( 2012l j called MOOD (multi-dimensional optimal order detection), and which adopts an 
a posteriori approach to the problem o f limiting of high orde r schemes in the finite volume framework. In a few words, the 
novel a posteriori DG limiter method of Dumbser et al. f (I 2014 I) consists of (i) computing the solution by means of an unlimited 
ADER-DG scheme, (ii) detecting a posteriori the troubled cells by applying a simple discrete maximum principle (DMP) and 
positivity of density and pressure on the discrete solution, (iii) creating a local sub-grid within these troubled DG cells, and 
(iv) reeomputing the discrete solution at the sub-grid level via a more robust Total Variation Diminishing (TVD) or Weighted 
Essentially Non Oscillatory (WENO) finite volume scheme. The final non-oscillatory DG solution on the main grid is then 
recovered from the subcell averages by means of a finite-volume reconstruction operator that acts on the cell averages of the 
subgrid. In the present paper we apply this idea for the first time to solve the RMHD equations in combination with space- 
time adaptiv e mesh refinement a nd time-accurate local time stepping, extending a sim ilar work proposed for classical fluid 
dynamics by 


Zanotti et al. 


2 OI 4 I: iFechter fc Miin7j|2ni5h . 


nentand time-accurate local time stepping, extending a sim ilar worK proposed tor classical nuid 
( 2015 ). For alternative work on DG subcell limiters see also ( Gasoni et al. boiffi Sonntag V Munzl 


The plan of the paper is the following. In Section [2] we report the RMHD equations and the basic physical assumptions, 
while Section [ 3 ] is devoted to the presentation of the numerical method, which is validated in Section |4l Section [5] contains a 
first simple analysis of the turbulence induced by the Kelvin-Helmholtz instability, while in Section [6] we conclude the work. 
We have considered a flat spacetime in pseudo-Gartesian coordinates, namely the metric = diag(—1,1,1,1), where Greek 
letters run from 0 to 3 and Latin letters i^j^ k,... run from 1 to 3. The speed of light is set to c = 1 and we make use of the 
Lorentz-Heaviside notation for the electromagnetic quantities, such that all factors disappear. Finally, we use Einstein 
summation convention over repeated indices. 
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2 MATHEMATICAL FORMULATION AND PHYSICAL ASSUMPTIONS 

The energy-momentum tensor of a single-component plasma with infinite conductivity is given by (lAnilelll99(][ l 

T'*" = {ph + E)u'^u‘' + {p + - bn-', (1) 

where is the four velocity of the fluid, is the four vector magnetic field, = b^b^, while h, p and p are the specific 
enthalpy, the rest mass density and the thermal pressure, each of them measured in the co-moving frame of the fluid. The 
metric of the spacetime is the Minkowski one, namely = diag(—1,1,1,1). We recall the in ideal MHD the electric 

field in the comoving frame of the fluid vanishes. If we instead select a static laboratory observer defined by the four-velocity 
vector = (— 1 , 0 , 0 , 0 ), then the electric field and measured in such a frame are related to the electromagnetic tensor 
and to its dual by 


F^^ = - E^n + 

F*^" = -B^n^ -e 






Bx'Gh 
E\ n^ , 


( 2 ) 

(3) 


where is the completely antisymmetric spacetime Levi-Civita tensor, with the convention that = 1. Note that the 

four vectors of the electric and of the magnetic field are purely spatial, i.e. E^ — B^ — 0, E^ — Ei^ B^ — Bi. Moreover, the fluid 
four velocity and the standard three velocity in the laboratory frame are related as v'^ — jW^ where W = {1 — 
is the Lorentz factor of the fluid. We stress that the electric field does not need to be evolved in time through the Maxwell 
equations, since within the ideal MHD assumption it can always be computed a posteriori as F = — T x F. In the rest of the 
paper we also assume that the fluid obeys the ideal gas equation of state, namely 


P = Pe(7 - 1): 


(4) 


where e is the specific internal energy, which is a function of the temperature only, and 7 is the adiabatic index. The equations 
of ideal RMHD, which in covariant form are 


= 0 , 

V„T“^ = 0, 

VaF*“^ = 0, 

for numerical purposes are better expressed in conservative form as ( Komissarov 19991 : Balsara 2001a) 

dtu + diE = 0 , 

where the conserved variables and the corresponding fluxes in the i direction are given b }0 


D 


FF 

Sj 

, r = 

wi 

u 


B^ 




(5) 

( 6 ) 

(7) 

( 8 ) 


(9) 


The conserved variables (F, Sj, U, B^) are related to the rest-mass density p, to the thermal pressure p, to the fluid velocity 
Vi and to the magnetic field B^ by 


D = pW, 

Si = phW^Vi + CijkEjBk, 

u = phw^ - p + + B^), 

where Cijk is the spatial Levi-Civita tensor and Sij is the Kronecker symbol. The spatial tensor W] in 
momentum flux density, is 


Wij = phW^ViVj — EiEj — BiBj + 


p+i{E^ + B^) 


( 10 ) 

(11) 

( 12 ) 

representing the 


(13) 


where Sij is the Kronecker delta. Eqs. © above include the divergence free condition V • F = 0. Although the Maxwell 
equations guarantee that such a constraint is mathematically fulfilled for all times if it is satisfied in the initial conditions, 
from a numerical point of view specific actions must be taken in order to preserve the divergence-free property of the ma gnetic 
field d uring the evolution of the system. Several strategies have been proposed over the years to solve this prob lem [see 


Toth 


(l200(]ll for a review]. In this paper we have adopted the so called divergence-cleaning approach presented in IPedner et al. 


^ Although f ormally written in conservativ e form, the evolution of the magnetic field is based on Stokes’ theorem rather than on Gauss’ 
theorem. See iLondrillo Sz Del Zann^ (l2000li for a careful discussion about these aspects. 
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(l2002h . which amounts to augmenting the system (|8]) with an additional equation for a scalar field in order to propagate 
away the deviations from V • B = 0. Hence, we must solve 

dt^ + diB^ = , (14) 

while the fluxes for the evolution of the magnetic field are also modified, namely —)► The damping 

coefficient k in Eq. d drives the solution towards V ‘ B — 0 over a tim e scale l/n. In our calculations we have typically used 


n E [1 ; 10]. More details about this approach can be found in Komissarov (2007), Palenzuela et al. ( 2009h . Dionysoponion et al 

(l20l4 . 

As well known, in the relativistic framework the conversion from the conserved variables (D, Si,U, Bi) to the primitive 
variables (p, p, Vi, Bi), which are needed for the computation of the fluxes, is not analytic, and a numerical root-finding ap proach 
is therefore needed. In our numerical code we adopted the third method reported in Sect. 3.2 o f Del Zanna et al.l (|2007l ) . A full 
account about alternative methods to invert the system (nnii-dni) was given i 


Noble et a] 


I] (l200dl 


about the mathem a tical properties o f the RMHD equatio n s can be found in Balsara fc Snicerl dlflflfll i: Komissarov ( IQQflh : 


Additional information 


Anton et al. ( 2006l l : Del Zanna et al. ( 2007 b Anton et al.l (|2010h . The latter, in particular, contains a detailed discussion 


about the renormalization of the eigenvectors of the associated Jacobian. 


3 NUMERICAL METHOD 


The numerical scheme that we adopt results from the combination of a few different steps, which in principle could be used 
separa tely. Here we provide a brief but self-consistent presentation of the scheme, while addressing to iDumbser et al.l (|2014r ) 


and to Zanotti et al 


p rovide 

I (|2015l l 


for additional discussion. 


3.1 Basic mathematical definitions 


We use spatial Cartesian coordinates over a domain Q which is composed by elements Ti as 

Ne 

n=\jTi, 

i=l 

where the index i ranges from 1 to the total number of elements Ne- At the generic time E 
is represented within each cell Ti by polynomials of maximum degree A ^ 0, namely 


(15) 


the numerical solution of Eq. 


Ufc(x,C) = = $i(x)u" X € Ti 


(16) 


1=0 


where u/, is referred to as the discrete representation of the solution, while the coefficients u" are the degrees of freedom^ 
In one spatial dimension, the basis functions ^i(x) are given by the Lagran ge interpolation polynomials, all of degree N, 
which pass through the (N + I) Gauss-Legendre quadrature points (|Solinll200^ ) . The resulting basis is therefore a nodal basis, 
with the property that ^i(xk) = Sik, where Xk are the coordinates of the Gauss-Legendre nodal points. In multiple space 
dimensions, the basis functions 4>z(x) are the dyadic products of the one-dimensional basis. 


3.2 The Discontinuous Galerkin scheme 


The system of equations (|8j) i s in conservative form and, therefore, numerical schemes derived from it are guaranteed to 
converge to the weak solution ( Lax fc Wendrofj I960l b even if this contains a discontinuity. The vast majority of numerical 
schemes for the solution of the RMHD equations use either conservative finite difference or finite volume schemes, which 
incorporate in a natural way Riemann solvers, thus assuring the upwi nd property of the method. An effective alternative 
to these approaches is represented by Discontinuous Galerkin methods ( Gockburn fc Shulll99ll . Il989l : iGockburn et al.l[l990l : 
Gockburn fc Shulll998l i. which still exploit the conservative form of the equations and the usage of Riemann solvers, but which 
evolve in time the degrees of freedom with respect to the chosen basis, rather then the point values or the cell averages of the 
solution. 

To illustrate the method, we first multiply the governing equations (|8]) by a test function identical to the spatial 


^ Here we slightly abuse of the Einstein summation convention, which is adopted even if the mute indices over which the summation is 
performed do not refer to co-variant and contra-variant vectors. 
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basis functions of Eq. iisi. After that, we integrate over the space-time control volume Ti x If we integrate by parts 

in space the flux divergence term, we obtain 

J j ^k^^dyidt+ J j {^h) - ^dSdt - J j • f (u/^) dxdt = 0, (17) 

Ti dTi Ti 

where n is the outward pointing unit normal vector on the surface dTi of element Ti. The second term of Eq. ca contains a 
surface integration, which is conveniently performed through the solution of a Riemann problem at the element boundary, like 
in traditional conservative finite volume schemes. This guarantees that the final method is upwind. The time integration 


of Eq. (pTl) can be performed through Runge-Ku t ta schemes, leading to RKDG schemes ( Cockburn fc Shu 1991 . 19891 : 


Cockburn et al. 198! 


1990 


Cockburn fc Shu 


19981: 


H.Zhu J.Qiu 2013l i but it can also be obtained through the ADER 


philosophy, see (jPumbser fc Mun'z 20061 : Qiu et all 2005l i. More precisely, we want to devise a one-step time integration 


scheme for ina, while preserving high order of accuracy both in space and in time. This can be done, provided an approximate 
predictor state is available at any intermediate time between and and with the same spatial accuracy of the initial 
DC polynomial. We have denoted this spacetime predictor solution c\h with a separate symbol, to distinguish it from the 
discrete solution of the DC scheme Uh- After inserting Uh, as given by ca, in the first term of (Ha and by using the spacetime 
predictor q/^ in the other terms, we find the following one-step ADER discontinuous Galerkin scheme: 

( \ 

J - uf) + J J ^kG {cih ,qt) ■ ndSdt - J J V'I’fe • f (q/i) dxdi = 0 . 


(18) 


In the equation above ^ is a numerical flux function, which in practice is given by a Riemann solver, while q^ and q^ are 
the corresponding left and right states of the spacetime predictor solution that is typical of ADER schemes. This will allow 
us to compute the spacetime integrals of the second and of the third terms of Eq. (flSl) to the desired order of accuracy. The 
strategy for obtaining the predictor q^ from the DG polynomials Uh is explained in SectjTS] Concerning the choice of the 
Riemann solver, in this paper we have used the simple Rusanov flux and the HLL solver dTorolIlQQflh . 


3.3 The spacetime discontinuous Galerkin predictor 

In the original ADER approach by Tit are v fc Torol ( 2002l i and Tit are v fc Torol ( 2005l i . the time evolution q^ of the data u^ 
available at time is obtained by means of the so-called Cauchy-Kowalevski procedure, which implies a Taylor expansion 
in time, and a subsequent replacement of time derivatives with spatial derivatives through the governing system of PDEs. 
As simple as it is in principle, this approach becomes prohibitively complex for highly no n-linear systems of equations. It 
has been successfully implemented for the classical Euler equations (Dinnbser^^^ j2007 i but it has never been extended 
to the relativistic regime. In the modern ADER version proposed by Dumbser et ali (|2008l b the time evolution is instead 
performed trough a spacetime discontinuous Galerkin predictor, which operates locally for each cell. To illustrate the method, 
we first transform the PDE system of Eq. (|8]) into a space-time reference coordinate system (^, 77 , C?Hence, the space-time 
control volume Cijkn — 1 ] x [yj_i;yj_^i] x x is mapped into the space-time reference element 

Te X [0; 1] as 


x = y = y.i TyAyj, z = Zj^_i (Azk, t = t"" + rAt, 


where Te = [0; 1]^ denotes the spatial reference element in d spatial dimensions. In these reference coordinates, Eq. 
rephrases into 

dr 

where 


+ ve-r (u) = o. 


(19) 

® 

( 20 ) 


f* := Ai ■ f(u) 


( 21 ) 


with ^ = {i,ri,Q) and Vj = d^/dx ■ V. We now multiply (I20II by a space-time test function 9k = dk{^,T) and integrate over 
the space-time reference control volume Te x [0; 1], to obtain 


■'«* + // 

0 Te 0 Te 


Ok'V^ • (u) d^dr = 0. 


( 22 ) 


The discrete spacetime solution of equation (|2^ is the q^h that we have mentioned above. In analogy to Eq. (IlGl) . we expand 
it as 


q/> = cih{^,T) = Oicii. 


( 23 ) 
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Something similar is done for the fluxes, which are represented as 

K = il{i,T) = eii:. (24) 

Both the space-time test function Ok in Eq. (|^ and the basis functions 0i are chosen as dyadic products of Lagrange 
interpolation polynomials passing through the Gauss-Legendre quadrature points. As a result, the degrees of freedom for the 
fluxes can be computed as the point-wise evaluation of the physical fluxes, i.e. 

fr=f*(qO- (25) 


Since the time evolution of the discrete solution is now hidden in the basis functions, we can integrate the first term by 
parts in time in (I22p . which allows us to introduce the DG solution as initial condition at time in a weak form. 

We thus obtain 


i i 

J J ek{^,0)uhd^- J J^qHd^dT + j j ft didT = Q. 

Te Te Te 0 Te 


(26) 


Inserting (|23l) and (|2^ into Eq. (l26l) provides ( Dumbser et al. 2008 : Hidalgo Dumbser 2011 : Dumbser Zanotti 2009l i 



which is a nonlinear system to be solved in the unknown expansion coefficients q^. A few comments should be given at this 
stage. The first one is that, being local in space, the spacetime discontinuous Galerkin predictor does not require the solution 
of any Riemann problem, which is instead invoked in the global scheme (HHi. The second comment is that the discontinuous 
Galerkin predictor just described can be used a lso in combination w i th mo re traditional finite volume schemes, which has been 
done for the RMHD equations for instance in Zanotti fc Dumbser (l2015lV The third comment is that this approach, unlike 
the original ADER approach, remains v alid even in the presence of s tiff source terms, as it has been done f or var ious physical 
systems by Dumbser Zanotti ( 20Q9l j : Hidalgo fc Dumbser (2011); Zanotti et al. (2011b Dumbser et al. (2012). Einally, we 
emphasize that for DG schemes the timestep must be restricted as (|Krivodonova R.Qin 2013h 


Ai < t 


1 


h 


where h and |An 


d(2iV + l) |A^ax| ’ 

are a characteristic mesh size and the maximum signal velocity, respectively. 


(28) 


3.4 An a posteriori subcell limiter 


Should we implement the ADER-DG scheme as it is described in the two previous Sections, we would obtain a numerical 
scheme capable of resolving smooth solutions with an order of accuracy equal to A + 1, where N is the degree of the 
chosen polynomials, but totally inadequate for discontinuous solutions, for which the Gibbs phenomenon would quickly lead 
to spurious oscillations and even t o a breakdown of the scheme. A novel idea for an a posteriori limiter has been recently 
proposed by Dumbser et al. I ll2014[l and it works as follows. 


• The unlimited ADER-DG scheme GUI is first used to evolve the solution from time to producing a so-called 

candidate solution inside each cell. 

• The candidate solution is then checked against two different criteria to verify its validity, namely 


(i) Physieal admissibility deteetion: if the conversion from conservative to primitive variables fails, or if either the pressure 
or the rest mass density drops below a threshold value, or if we encounter superluminal velocities, then the cell is flagged 
as troubled. 

(ii) Numerieal admissibility deteetion: if the polynomial representing the candidate solution does not lie between the mini¬ 
mum and the maximum of the polynomials representing the solution at the previous time step in the set Vi, then the cell 
is flagged as troubled. The set Vi contains the cell Ti and all its Voronoi neighbor cells that share a common node with Ti. 
This second detection criterion is specifically designed to remove spurious Gibbs oscillations from the solution. 

• As soon as a cell is flagged as troubled at the future time it generates a local sub-grid formed by Ns — 2N + 1 cells 
per space dimension, each of which is assigned a subcell average v/i(x,t^) by means of a L 2 projection obtained from the DG 
polynomial at the previous time level i.e. 

= id—I / Uh(x,i’")dx = 1-2—- / ut4>i{x)dx, VS'ij € 5,, (29) 

\^Fj\ JSij J Sij 

where Si = |Jj ^^e sub-grid cells. In this way the high accuracy of the DG polynomial is transferred to the 

subgrid level before the spurious oscillations arise. We have chosen Ns = 2N + 1 in order to guarantee that the maximum 
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timestep of the ADER-DG scheme on the main grid (c.f. Eq. (l28l) ) matches the maximum possible time step of the ADER 
finite volume scheme on the sub-grid. 

• The alternative data representation, provided by Eq. (123), is now used as initial condition to evolve the discrete solution 
with a more robust finite volume scheme on the sub-grid. This is done by resorting to either an ADER-WENO finite volume 
scheme, or to an even more robust second order TVD shock captu r ing scheme. Eor detai l s abou t the implementation of WENO 
within our ADER framework we refer to Dumbser et al. ( 2Q13l i: Zanotti fc Dumbser ( 2015h . In practice, on the sub-grid a 
new evolution from time to is performed combining a third order WENO finite volume scheme with the spacetime 
discontinuous Galerkin predictor described in Sect. 13.31 Only for particularly challenging problems we sacrifice WENO in 
favor of a simpler second order TVD scheme. We emphasize that both the DG scheme on the main grid as well as the WENO 
finite volume scheme on the sub-grid are one-step ADER schemes. 

• The last step requires that the new solution at time over the sub-grid is projected back to the main grid. This is 
done imposing that 


J u/i(x, = J v/,(x,t^^^)dx, WSij 


G Si 


(30) 


which is a standard reconstruction problem in hi gh order finite volume methods ( Barth fc EredericksonI 199C : Tit arev Torol 


2 OO 4 I : Titarev fc Torol 20051 : Dumbser et al. 2013l i and spectral finite volume schemes 


In all the numerical simulations described in Sect.lH the DG scheme over the main grid has been implemented with N G [2; 5], 
hence up to the sixth order of accuracy both in space and in time, while the WENO scheme on the sub-grid is always at the 
third order. Moreover, in the Eigs. [T] and [4| below we have represented in blue the unlimited cells, namely those that have 
been successfully evolved through the standard ADER-DG scheme, while we have represented in red the troubled cells, which 
required the activation of the subcell limiter. 


3.5 Adaptive Mesh Refinement 


The whole scheme described s o far can be implemen t ed over adapti v ely re fined meshes (AMR), together w ith time-accurate 
local time-stepping (LTS). In Dumbser et al. ( 2013l j. Zanotti et al. ( 2015l j and Zanotti V Dumbser (2015) we have already 
described all details of our AMR strategy, hence, in the following we rec all the most imp ortant aspects only. 

Our AMR approach can be referred to as a ’’cell-by-cell” refinement ( Khokhlov 1993), according to which every cell Ti is 
individually refined with no creation of grid patches. As customary, the refinement criterion involves up to the sec ond order 
derivative of a suitable indieator funetion in terms of which a suitable refinement function is build (|Lhnerlll987l i. 




i {d‘^^/dxkdxif 


\ Efc.i [ ( I d^/dxk li+i + I d^/dxk li ) / Ax, 


l +e 


^2 

dxkdxi 


14-1 


(31) 


The refinement function Xm is checked for each cell Ti, and if Xi > Xref, the cell is refined, while it is recoarsened if Xi < Xrec- 
In most of our simulations, except for Sect. ITT] where we have used ^ = By, the indicator function has been assumed to be 
the relativistic mass density, i.e. ^ = D = Wp. The implementation of the AMR infrastructure is based on the following 
general rules: 


• A maximum level of refinement imax is chosen, such that 0 ^ ^ ^ ^max, where i indicates the actual refinement level. 

• When a mother eell Ti is refined, it generates ehildren cells, where d is the spatial dimension, while r is the refinement 
factor, typically chosen between 2 and 4. 

• Each cell Ti, at any level of refinement, is given a specific status, denoted by a for convenience, with the following meaning 

(i) active cell ((7 = 0), updated through the standard ADER-DG scheme; 

(ii) virtual child cell {<7 — 1 ), updated according to standard L 2 projection of the high order polynomial of the mother cell 
at the {i. — l)-th level; 

(hi) virtual mother cell (a = —1), updated by recursively averaging over all children cells from higher refinement levels. 

A virtual child cell has always an active mother cell, while a virtual mother cell has children with status (7^0. 

• Only active cells (a = 0) can be refined. Hence, if a virtual cell needs to be refined, it must be first activated. 

• The levels of refinement of two cells that are Voronoi neighbor^ of each other can only differ by at most unity. Moreover, 
every cell has Voronoi neighbors, which can be either active or virtual, at the same level of refinement. 


3 


The Voronoi neighbors Vi of a cell Ti are cells which share common nodes. 
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2D circularly polarized Alfven Wave problem — ADER-DG-Pat + WEN03 SCL 



N, 

Li error 

L 2 error 

Loo error 

Li order 

L 2 order 

Loo order 

Theor. 

(N 

30 

2.9861E-3 

7.2314E-4 

4.2388E-4 

— 

— 

— 


1 

60 

2.9229E-4 

8.0346E-5 

7.0230E-5 

3.35 

3.17 

2.59 

Q 

Q 

90 

8.8069E-5 

2.5059E-5 

2.3319E-5 

2.95 

2.87 

2.72 



120 

3.6687E-5 

1.0900E-5 

1.0948E-5 

3.04 

2.89 

2.63 


CO 

15 

1.2671E-4 

2.5939E-5 

1.1433E-5 

— 

— 

— 


1 

20 

3.1455E-5 

6.5949E-6 

2.9456E-6 

4.48 

4.76 

4.71 

A 

Q 

25 

1.1743E-5 

2.5410E-6 

1.4527E-6 

4.41 

4.27 

3.17 

4 


30 

5.7046E-6 

1.2767E-6 

7.5875E-7 

3.96 

3.77 

3.56 



10 

6.6600E-5 

1.4648E-5 

7.5420E-6 

— 

— 

— 


1 

15 

7.8640E-6 

1.9384E-6 

1.2828E-6 

5.26 

4.98 

4.36 

IT 

Q 

20 

1.8748E-6 

4.9562E-7 

3.6520E-7 

4.98 

4.74 

4.36 

0 


25 

6.1631E-7 

1.6408E-7 

1.3283E-7 

4.98 

4.95 

4.53 



Table 1. Li,L 2 and Loo errors and convergence rates for the 2D circularly polarized Alfven wave problem for the ADER-DG-Piv 
scheme with subcell limiter and adaptive mesh refinement. Two levels of refinement have been used with a refinement factor r = 3. The 
errors have been computed for the variable . 


The rules above are quite general, and they would still hold even if a pure finite volume scheme was adopted. In addition to 
them, a few more instructions are needed when the AMR framework is combined with the presence of the limiter for the DG 
scheme. Namely, 

• The virtual children cells inherit the limiter status of their active mother cell. 

• If at least one active child is flagged as troubled, then the (virtual) mother is also flagged as troubled. 

• Cells which need the subcell limiting cannot be recoarsened. 


A proper description of th e AMR-pro.iection a nd of the AMR-averaging at the sub-grid level involving different levels of 


refinement can be found in I Zanotti et al.l (|2015l l . 


4 NUMERICAL TESTS 


4.1 Convergence test 


We have tested the convergence of our new nume rical scheme by considering the propagati on of a circularly polarized Alfven 


wave, for which an analytic solution is known (|Komissarovl 119971 : iDel Zanna et al.l 120071 ) . Choosing x as the direction of 


propagation, and 77 as the amplitude of the wave, the magnetic field is given by 


Bx 

= Bo 

(32) 

By 

= pBo cos[A;(x — VAt)] 

(33) 

Bz 

= pBo sin[/c(x - VAt )], 

(34) 


where Bq is th e uniform magne t ic fiel d along x, /c is the wave number, while va is the Alfven speed at which the wave 
propagates fsee lPel Zanna et al. ( 2007 ) for its analytic form). The vector tips of the transverse velocity field describe circles 
in the yz plane normal to Ro, according to 

Vy = -VAByjBQ, Vz = -vaBzIBo . (35) 


We have used p = p = Bo = r]=l, and since the wave is incompressible, the background values of p and p are not affected. 
The test has been performed in two spatial dimensions, using periodic boundary conditions, over the computational domain 
LL = [0; 27r] X [0;27r]. We compare the numerical solution with the analytic one after one period T = LjvA = 27tIva- The 
results of this analysis are reported in Tab.[T] which report the Li, L 2 and Loo norms of the error of B^. The Rusanov flux has 
been adopted, with ^max = 2 and a Courant factor CFL = 0.8. We emphasize that, due to the smoothness of the solution, the 
subcell limiter is never activated. As it is apparent from the table, the nominal order of convergence is essentially confirmed. 


4.2 Riemann problems 


Once the convergence properties have been verified, we consider a few relevant shock-tube problems to test the n ew ADER- 
DG-AMR method. Specifically, we concentrate on two classical Riemann problems for RMHD, already proposed bv IvanPutten 
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Figure 1. 3D view of the density variable and the corresponding AMR grid. Top panel: RPl at tfinal = 0-4 (coarsest grid of 40 x 5 
elements). Bottom panel: RP2 at tfinal = 0-55 (coarsest grid of 25 x 5 elements). The limited cells, using the subcell ADER-WEN03 
finite volume scheme, are highlighted in red, while unlimited DG-P 3 cells are highlighted in blue. 
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Problem 


P 

{Fx 

Vy 

Vz) 

P 

{B^ 

By 

Bz) 

^final 

7 

RPl 

a: > 0 

0.125 

0.0 

0.0 

0.0 

0.1 

0.5 

-1.0 

0.0 

0.4 

2.0 

iTest 1 in Balsara r2001aB 

a: ^ 0 

1.0 

0.0 

0.0 

0.0 

1.0 

0.5 

1.0 

0.0 

RP2 

a: > 0 

1.0 

-0.45 

-0.2 

0.2 

1.0 

2.0 

-0.7 

0.5 

0.55 

5/3 

(Test 5 in Balsara f2001aB 

a: ^ 0 

1.08 

0.4 

0.3 

0.2 

0.95 

2.0 

0.3 

0.3 


Table 2. Initial conditions for the one-dimensional Riemann problems. 






Figure 2. RPl: physical variables interpolated along a ID cut on 200 equidistant points at tfinal = 0-4, starting from a coarsest grid of 
40 X 5 elements by using the ADER-DG-P 3 scheme supplemented with the a posteriori ADER-TVD subcell limiter. 


(|l993h and classified as Test 1 and Test 5 in Table 1 of Balsaral (2001a). The corresponding initial conditions, referred to as 


RPl and RP2 in the following, are given in Tab.O which reports also the final times and the adiabatic indices|f| 

The two chosen Riemann problems are solved along two coarse grids of 40 x 5 and 25 x 5 elements, respectively. Then, 
the intial grid is adaptively refined in space and time according to r = 3 and ^max = 2. The computational domain is only 
formally two-dimensional, since the second direction y acts as a passive one. Both tests have been solved using the ADER- 


^ Th e adiabatic index 7 of RPl in unphysical, as it violates Tau b’s inequality base d on kinetic theory (|Taublll94^ : lMignone fc McKinnevI 
l20Q7h but it is fixed equal to 2 anyway to ease comparison with IvanPnttenI (Il993h and iBalsaral (l 2001 alh 
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Figure 3. RP2: physical variables interpolated along a ID cut on 200 equidistant points at tfinal = 0.55, starting from a coarsest grid 
of 25 X 5 elements by using the ADER-DG-P 3 scheme supplemented with the a posteriori ADER-WEN03 subcell limiter. 


DG-P 3 scheme, but they differ in the subcell limiter, which is the second order TVD finite volume scheme for RPl, while it 
is the third order ADER-WENO finite volume scheme for RP2. The HLL solver has been used for both RPl and RP2 with 
a Courant factor CFL = 0.5. The damping factor for the divergence-cleaning procedure is set to /c = 10. 

Figure [1] shows the three-dimensional plot of the solution for the rest mass density and the corresponding AMR grid, by 
plotting the real DG polynomials (highlighted in blue) for every single unlimited cell and the piecewise linear interpolation of 
the ADER-WENO limiter along the subcell averages (highlighted in red) fo r the limited cells. A reference solution for these 
Riemann problems is computed with the exact Riemann solver proposed by Giacomazzo fc Rezzolla ( 2006h . Eigures[2]and[3] 
show the comparison with the reference solution by plotting the rest mass density, the x— and the y— velocity components 
and the y— component of the magnetic field, interpolated over a one-dimensional cut composed of 200 equidistant points at 
the final state. A remarkable agreement between the numerical and the reference solution is obtained. All the waves are well 
captured, five for RPl and seven for RP2. More specifically, RPl has a left-going and a right-going fast rarefaction wave, a 
left-going compound wave, a central contact discontinuity, and a right-going slow shock. RP2 has instead a left-going and a 
right-going fast shock, a left-going and a right going Alfen wave, a left-going rarefaction wave, a central contact discontinuity 
and a right-going slow shock. Due to the combined action of the subcell-limiter and of AMR, all discontinuities are resolved 
within just one cell or two cells at most. We note that, while the compound wave is absent by construction in the exact 
solution, its width in the numerical solution is rather small and its amplitude is also comparatively smaller with respect to 
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that obtaine d with other numerica l schemes, indicating that this might really be a numerical artifact. However, see also the 

(l2009ll 


discussion iniMignone et al 


The small asymmetries visible in Fig. [T]for RPl along the passive y direction are attributable to the joint interaction 
between: (1) the lack of reconst ruction in cha racteristic variables, which could typically help in these cases; (2) some residual 
post-shock oscillations, that in ( Balsara IQQSl l were suppressed by means of artificial viscosity. In spite of these small defects, 
these results show the capabilities of the new scheme, which does not resort to any artificial viscosity, in resolving the strongly 
non-linear waves of RMHD equations, for which an unlimited DG schemes would catastrophically fail. 


4.3 The rotor problem 

As a first genuinely two d i mensi onal test we consider the relativistic version of the MHD rotor problem, origina l ly pro ¬ 
posed by Balsar a &: Spicer ([l999 ), and so l ved b y a n umber of autho r s over the years, including Del Zanna et al. (2003), 
Dumbser Zanottil (|2009h . iLoubere et al. (|2014h and Kim Balsaral (|2014l i. The computational domain is chosen to be 


Q = [—0.6, 0.6] X [—0.6, 0.6], discretized on a coarse initial grid formed by 40 x 40 elements. The AMR framework is activated 
with a refinement factor r = 3 and a number of refinement levels ^max = 2. In this problem a cylinder of a high density fluid is 
rotating rapidly with angular velocity cu, surrounded by a low density fluid at rest. The initial conditions are in fact given by 

1.0 

1 ^ 0 P = (36) 


which imply an initial maximum Lorentz factor VFmax ~ 10 at r = 0.1. Transmissive boundary conditions are applied at 
the borders. The spinning of the rotor produces torsional Alfven waves that are launched outside the cylinder, transferring 
amounts of its initial angular momentum into the external medium. The simulation is performed without any linear taper, 
that means the physical variables between the internal rotor and the fluid at rest are really discontinuous. The adiabatic index 
is 7 = 4/3. For this test, the P5 version of our ADER-DG scheme was used, combined with the Rusanov Riemann solver. Due 
to the challenging nature of the problem, a robust second-order TVD scheme, rather then the standard WENO scheme, has 
been used on the subgrid where the limiter is activated. 

Eig.[4|shows the rest-mass density, the thermal pressure, the relativistic Mach number M and the magnetic pressure pMag 
at time t = 0.4. The latter are computed according to 


M = 


Wv 


_ 17 2 _ 

PMag ~ 2^ ~ 


B‘^/W + (v • B) 


WsVs' 2 " 2 ’ 

where Vg is the speed of sound and Wg = (1 — is the corresponding Lorentz factor. Although an analytic solution 

is not available for this test, the results shown are in very good qualitative agreement with those already reported in the 
literature. In particular, the maximum Lorentz factor of the rotor, which is considerably slowed down by magnetic braking, is 
Wmax ~ 2.1. Moreover, the adopted divergence-cleaning approach works accurately as expected, with no appreciable spurious 
oscillations generated in the rest mass density or in the magnetic field. Lastly, the behaviour of the space-time AMR and of 
the a posteriori limiter is depicted in the two bottom panels of Eig. [D the final mesh is shown in the left panel, whereas in 
the right the troubled zones are represented in red. Glearly, the activation of the limiter becomes necessary only in a limited 
number of cells, and precisely where discontinuities are stronger. 


4.4 Cylindrical blast wave 


As a second two dimensional academic test we have considered the cylindrical expansion of a blas t wave in a plasma with 


an initially uniform magnetic field. This is not oriously a severe test, which became canonical a fter 


it has been solved by several authors, including iLeismann et al 


Del Zanna et al.l (|2007l i ; iDumbser fc Zanottil ([21 


Komissaro^ 



and 


The initial conditions are prescribed by assuming that, within a radius R = 1.0, the rest-mass dens ity and the pressure are 


p — 0 .01 and p = 1, while outside the cylinder p = 10 ^ and p = 5 x 10 Like in Komissarov ( IQQOl i and in Del Zanna et al, 


( 20071 1. the inner and outer values are joined through a smooth ramp function between r = 0.8 and r = 1, to avoid a sharp 
discontinuity in the initial conditions. The plasma is initially at rest and subject to a constant magnetic field along the 
x-direction. In our tests we h ave consider e d two different magnetizations, the first one with Bx — 0.1, corresponding to the 
intermediate value chosen bv iKomissaro 3 lll999lV and the second one with Bx = 0.5. We have solved this problem over the 
computational domain Q — [—6,6] X [-6, 6], with 40 X 40 elements on the coarsest refinement level, r = 3 and ^max = 2. We 
have used the Rusanov Riemann solver with the P 3 version of the ADER-DG scheme. Also for this test, a robust second-order 
TVD scheme has been used on the subgrid where the limiter is activated. The results for Bx — 0.1 are shown in Eig. O 
which reports the rest-mass density, the thermal pressure, the Lorentz factor and the magnetic pressure at time t = 4.0. The 
wavepattern of the configuration at this time is composed by two main waves, an external fast shock and a reverse shock, 
the former being almost circular, the latter being somewhat elliptic. The magnetic field is essentially confined between them. 
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Figure 4. Solution of the RMHD rotor problem at time t = 0.4, obtained with the ADER-DG P 5 scheme supplemented with the a 
posteriori second order TVD subcell limiter. Top panels: rest-mass density (left) and thermal pressure (right). Central panels: Mach 
number (left) and magnetic pressure (right). Bottom panels: AMR grid (left) and limiter map (right) with troubled cells marked in red 
and regular unlimited cells marked in blue. 
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Figure 5. Solution of the RMHD blast wave with Bx = 0.1 at time t = 4.0, obtained with the ADER-DG P 3 scheme supplemented 
with the a posteriori second order TVD subcell limiter. Top panels: rest-mass density (left) and thermal pressure (right). Central panels: 
Lorentz factor (left) and magnetic pressure (right), with magnetic field lines reported. Bottom panels: AMR grid (left) and limiter map 
(right) with troubled cells marked in red and regular unlimited cells marked in blue. 
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Figure 6. Solution of the RMHD blast wave with Bx = 0.5 at time t = 4.0, obtained with the ADER-DG P 3 scheme supplemented 
with the a posteriori second order TVD subcell limiter. Top panels: rest-mass density (left) and thermal pressure (right). Central panels: 
Lorentz factor (left) and magnetic pressure (right), with magnetic field lines reported. Bottom panels: AMR grid (left) and limiter map 
(right) with troubled cells marked in red and regular unlimited cells marked in blue. 
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while the inner region is almost devoid of magnetization. We have detected a maximum Lorentz factor Wmax ~ 4.3 along the 
X axis just on the back of the reversed shock. The two bottom panels show the AMR grid and the map of the limiter, which is 
activated along the two main shock fronts. In Fig. [ 6 l on the other hand, we have reported the results obtained for Bx = 0.5, 
again at t = 4.0. In this case, the external circular fast shock, which is visible in the rest-mass density and in the magnetic 
pressure, is very weak, while the magnetic confinement of the plasma is increased. The maximum Lorentz factor detected in 
this case is VFmax ~ 2 . 8 . 


4.5 Orszag-Tang vortex system 


Next, we have chosen the rel ativistic version of the wel l kno wn Orszag-Tang vortex p roblem, proposed by Orszag Tand 
(Il979l) . and later considered bv IPicone Sz Dahlburgl dlQOlh an d lPahlburg Picon3 (jlQSOl L The resistive case of this relativistic 


MHD problem has been investigated bv iDumbser Zanottil (|2009l L The initial conditions are given by 


1 3 . 

1 ,-^ sm 

4^/2 

while the adiabatic index is 7 = 4/3. The equations are discretized over the computational domain Q = [0, 27r] x [0,27r], 
with 30 X 30 elements on the coarsest refinement level at the initial state. Periodic boundary conditions are imposed at the 
borders and the Rusanov Riemann solver is adopted. The relevant AMR parameters are r = 3 and ^max = 2. We note that 
the maximally refined AMR mesh corresponds to a uniform grid formed of 270 x 270 = 72, 900 elements. Moreover, the P 5 
version of the ADER-DG scheme that we have adopted uses 6 degrees of freedom per spatial dimension, amounting to a total 
resolution of 2, 624, 400 spatial degrees of freedom. The computed solution for the rest-mass density is shown in the central 
column of Fig. [3 at times t = 0.5, 2.0, 3.0, 4.0 respectively. For comparison, the panels on the right column show the results 
of a simulation performed over the maximally refined uniform mesh, which can be used as a reference solution. Clearly, an 
excellent agreement between the AMR results and this reference solution is obtained. As before, this test confirms the ability 
of the proposed method for solving complex two dimensional problems and, by showing the critical cells which required the 
activation of the limiter, it provides an immediate visual sketch of the most delicate regions over the computational domain. 


(y) , sin (x), 0 , 1 , - sin (y), sin (2x), 0 


(38) 




5 THE RMHD KELVIN-HELMHOLTZ INSTABILITY 


A two-dimensi onal test that is not on ly academic but may be relevant to explain the observed phenomenology of extended 
radio-jets [see Martf fc Muller ( 20Q3l j and references therein], we consider the Kelvin-Helmholtz (K H) instability with an 
i nitial ly uniform magnetic field. Following the works of Mignone et al. ( 2009h . Beckwith fc Stone ( 201 ih and iRadice Rezzolla 
(| 2012 h . we choose the initial conditions as 


Vs tanh [{y - 0.5)/a] y > 0, 

-Vs tanh [{y + 0.5)/a] y ^ 0 , 


(39) 


where Vs = 0.5 is the velocity of the shear layer and a = 0.01 is its characteristic size. Although not necessary in principle, it 
is convenient to introduce a small transverse velocity to trigger the instability, hence fixing 

yoVs sin (27rx) exp [-{y - 0.5)^/cr] 2 / > 0 , 

-ryo^s sin (27rx) exp [-(y + 0.5)^/cr] y ^ 0, 


where r]o = 0.1 and a = 0.1. Finally, the rest-mass density is chosen as 

po + pi tanh [(y - 0.5)/a] 


P = 


po - pi tanh [(y + 0.5)/a] 


y > 0, 
y ^ 0, 


(41) 


with po = 0.505 and pi = 0.495. The adiabatic index is 7 = 4/3, the pressure is p = 1 everywhere, and we add a weak uniform 
magnetic field along the x— direction, namely Bx = 0.001. The simulations are run with the ADER-DG-P 3 scheme over the 
computational domain Q = [—0.5, 0.5] x [—1,1], using 50 x 100 elements on the coarsest refinement level at the initial state. 
Periodic boundary conditions are imposed along all borders and the Rusanov Riemann solver is adopted. AMR is activated 
with r = 3 and ^max = 2. In this simulation the solution on the subgrid has been evolved through a second order TVD scheme, 
which turned out to be more robust than the usual third order WENO method. Fig. [ 8 ] shows the rest-mass density field at 
various times, up to t = 30, and the corresponding development of the KH instability. Since no physical viscosity or resistivity 
is present, it is very difficult to judge about the physical nature of the tiny structures, especially secondary instabilities, 
which are produced during the evolutio n, and which have been shown to depend sensib l y on the order of accurac y of the 
scheme and on the Riemann solver used ( Beckwith fc Stone 2011 : Radice fc Rezzollal 2012 : Zanotti Dumbser 2015h . As the 
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Figure 7. RMHD Orszag-Tang vortex problem at times t = 0.5, t = 2.0, t = 3.0, t = 4.0, from top to bottom, obtained through the 
ADER-DG-P 5 scheme supplemented with the third order a posteriori ADER-WENO subcell limiter. Left panels: AMR-grid, troubled 
cells (red) and unlimited cells (blue). Central panels: Ps-solution obtained on the AMR grid. Right panels: Ps-solution obtained on the 
fine uniform grid corresponding to the finest AMR grid level. 
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Figure 8. RMHD Kelvin-Helmholtz instability at times t = 5.0, t = 10.0, t = 20.0, t = 30.0 from left to right, obtained through the 
ADER-DG-Ps scheme supplemented with the second order a posteriori ADER-TVD subcell limiter. The computed solution of density 
(top), AMR grid (center) and limiter map (bottom) are shown. 
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Figure 9. Power spectra for the RMHD Kelvin-Helmholtz instability at time t = 30.0 obtained through the ADER-DG-P 3 scheme. 


instability proceeds, the transition to a turbulent state occurs. Although our final time is not large enough to allow for a 
fully developed turbulent state, and a l thoug h this paper is not devoted t o a detailed study of relati vistic MHD turbulence 
(see instead the works by Zhang et al. ( 2009h : Zrake MacFadveu (2012); Garrison Nguv^ (2015)), we have nevertheless 
computed the power spectra of a few relevant quantities to confirm that the transition to turbulence is indeed taking place. 
Fig.lH in particular, shows the power spectra of the velocity field, of the pressure field and of the magnetic field, which have 
been computed according to 


Pvik) = 11 
^ J\k 


\k\=k 


|'i)(fc)|^ dk, 


Ppik) = 


L 


\k\=k 


\fi{k)\^dk, 


PB(fc) = 


/ 

J\h 


\k\=k 


\B{k)\^dk, 


(42) 


where k is the wave-number, while v{k)^ B(k), p{k) are the two-dimensional Fourier-transforms of v, p and B, respectively. 
For k ~ [20, 70], in the so-called inertial range where the dynamics of the turbulence is not affected by large sca le energy input s 
nor by dissipation, we approximately recover Kolmogorov’s trends, namely Pv{k) oc and Pp{k) oc k~^^^ ( Biskan ml2008 i. 

The power spectrum of the magnetic field is instead in qualitative agreement with results obtained by Izhang et al.l (|2009l i for 
fully turbulent configurations. A dedicated analysis to astrophysical RMHD turbulence will be presented in a separate work. 


6 DISCUSSION AND CONCLUSIONS 


We have proposed a novel approach for the numerical solution of the special relativistic magnetohydrodynamics equations, 
which is based on the discontinuous Galerkin finite element method, but with some crucial modifications. Pure DG schemes, in 
fact, cannot avoid the appearance of oscillations when discontinuities form in the solution. The common practice in these cases 
is to resort to either artificial viscosity, filtering, or to an a priori finite-volume-type limiting of the higher order moments of 
the DG scheme, thus spo iling the sub cell resolution capabilities of the DG scheme. On the contrary, following the recent works 
by Dumbser et al. ( 2014h and Zanotti et al. ( 2015h . it is possible to verify a posteriori the validity of the candidate solution 
provided by the pure DG scheme by applying a relaxed form of the discrete maximum principle and by checking the numerical 
solution for physical validity, i.e. for positivity, subluminal velocities and for possible failures in the the conversion from the 
conservative to the primitive variables. For those cells that violate any of the two criteria, we scatter the DG polynomials, 
computed at the previous (safe) time step, onto a set of Ng = 2N + 1 subcells along each spatial direction. After that, a 
traditional and more robust WENO finite volume scheme, or an even more robust TVD scheme, is applied at the subeell level, 
thus reeomputing the solution in the troubled cells. In this way, our special a posteriori subcell limiter of the DG method 
is intrinsically based on the governing partial differential equations, while standard DG limiters act independently of the 
governing PDF. The new, and thereby safe, subcell averages are subsequently gathered back into high order cell-centered DG 
polynomials on the main grid through a subgrid reconstruction operator. The a posteriori limiter that we have used can also be 
regarded as a discontinuity detector. This is clear, for instance, after looking at the bottom right panel of Fig. |4] and of Fig. (5] 
where the activation of the limiter occurs where strong discontinuities are present. In this respect, our a posteriori correction 
may resemble the logic behind a priori limiters based on shoek deteetors, several forms of which have been introduced, both 
in the Newtonian and in the relativistic framework. However, the a posteriori approach that we have developed has at least 
two evident advantages over traditional shoek deteetors. The first one is that it is both very simple and it applies unmodified 
for general systems of equations, while shock detectors become more and more elaborate as the complexity of the equations 
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increases^ The second advantage is that the a posteriori approach can capture all kind of discontinuous solutions, including 
contact discontinuities, which are invisible to shock detectors but often represent a serious challenge to the numerical scheme. 

We also stress that both the DG scheme on the main grid and the WENO (or TVD) finite volum e scheme on the subgrid 
are implemented with the local spacetime discontinuous Galerkin predictor of Dumbser et al. (|2008l b thus providing a high 


order one-step ADER scheme in time, with no need for the Runge-Kutta time discretization that is typically used in the so- 
called method of lines. Einally, adaptive mesh refinement (AMR) is used, implying that the two AMR operations of projection 
and averaging need to involve the subcell averages of the solution on the sub-grids. Eor particularly challenging tests and 
applications, we have found convenient to resort to traditional second-order TVD schemes on the subgrid, thus achieving, at 
least formally, the same robustness of those schemes. 

The new ADER-DG-AMR scheme, which has been validated over stringent academic tests, can contribute significantly to 
the numerical modeling of high energy astrophysics systems, such as extragalactic jets, gamma-ray bursts and magnetospheres 
of neutron stars. Work is already in progress to extend this approach to the full general relativistic regime. A particularly 
attracting field of application for the new method is represented by the study of relativistic turbulence, that we have considered 
here in a simplified and preliminary two-dimensional set up in which the transition to turbulence is induced by the Kelvin- 
Helmholtz instability. We plan to investigate this problem in the future by means of three-dimensional simulations in which 
the high order capabilities of DG schemes are fully exploited. 
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